Visualization of the template files associated with fix bond/react¶
https://docs.lammps.org/fix_bond_react.html
This notebook walks the user step-by-step to build a reactive simulation with a non-reactive force field using the fix bond/react command. Ethyl cellulose is used as an example.
In [1]:
from IPython.display import display, HTML
display(HTML("<style>.container { width:90% !important; }</style>"))
from polymaps.util import viz_lmp_template, build_react_template
from polymaps.util import read_lmp_data, write_lmp_data, viz_mol, viz_lmp_data
The two monomers before reacting¶
The LAMMPS data file can be generated with OPLS-AA force field from LigParGen: http://zarbi.chem.yale.edu/ligpargen
In [2]:
mass_df, df, bond_df, angle_df, dihedral_df, improper_df = viz_lmp_data('before.data', annotation=True)
pick the atoms to eliminate during polymerization: [28, 37, 53], new bond forms between [10, 38]¶
The selected atoms to be included in the template file should cover all atoms that are involved in the same 4-body, 3-body, 2-body interaction types defined by the force field.
all heavy atoms that are within the 3rd degree neighborhood of atoms [10, 37] and nearby hydrogen atoms should be selected
first we will build a adjacency list from the bond list, then a function tallies all neighbors within the 3rd degree connection.
In [3]:
adj_dict = dict(zip(df["id"].to_list(), [list(set(bond_df[bond_df["at1"]==x]["at2"].tolist() + \
bond_df[bond_df["at2"]==x]["at1"].tolist())) for x in df["id"]]))
initiator_atom_ids = [38, 9]
def find_neighbors_within_Nth_degree(N, atom_id, adj_dict):
neigh_list = []
if N < 0:
print("Error: N = " + str(N) + " < 0!!\n")
return []
curr_centers = [atom_id]
neigh_list = neigh_list + curr_centers
if N == 0:
return neigh_list
for currN in range(1, N+1):
curr_perimeter = []
for center in curr_centers:
new_neighbors = adj_dict[center]
curr_perimeter = list(set(curr_perimeter + list(set(new_neighbors) - set(neigh_list))))
curr_centers = curr_perimeter[:]
neigh_list = list(set(neigh_list + curr_perimeter[:]))
edge_atoms = []
for center in curr_centers:
new_neigh = adj_dict[center]
num_of_neighs_of_new_neigh = len(list(set(new_neigh) - set(neigh_list)))
if num_of_neighs_of_new_neigh > 0:
edge_atoms.append(center)
return neigh_list, edge_atoms
The selected atoms are determined separately, and then the atom ids are combined. Using the module function to build the reaction template file "rxn1_ec_unreacted.data_template"¶
In [4]:
template_name_before_reaction = "rxn1_ec_unreacted.data_template"
selected_atoms_mol1, edge_atoms_mol1 = find_neighbors_within_Nth_degree(3, 10, adj_dict)
selected_atoms_mol2, edge_atoms_mol2 = find_neighbors_within_Nth_degree(3, 37, adj_dict)
selected_atoms = sorted(list(set(selected_atoms_mol1 + selected_atoms_mol2)), reverse=True)
atom_id_remap = build_react_template(df, bond_df, angle_df, dihedral_df, improper_df, selected_atoms, template_name_before_reaction)
#viz_lmp_template(template_name_before_reaction)
template_df = df.loc[selected_atoms, :].copy(deep=True)
template_df["template_atom_id"] = template_df["id"].map(atom_id_remap)
Written to file: rxn1_ec_unreacted.data_template
initiator atoms (https://docs.lammps.org/fix_bond_react.html) are defined as atoms [2, 9] from 2 distinct molecules.
In [5]:
from polymaps.util import build_map_file
edge_atom_ids = edge_atoms_mol1 + edge_atoms_mol2
map_file_name = "rxn1_ec_map"
build_map_file(map_file_name,
template_df,
initiator_atom_ids,
edge_atom_ids,
[atom_id_remap[x] for x in [28, 37, 53]],
[atom_id_remap[x] for x in edge_atom_ids])
reaction template before the reaction¶
In [6]:
viz_lmp_template(template_name_before_reaction)
The two monomers after reacting¶
In [7]:
mass_df, df, bond_df, angle_df, dihedral_df, improper_df = viz_lmp_data('after1.data', annotation=True, double_bonds=[])
building template after the reaction¶
In [8]:
adj_dict = dict(zip(df["id"].to_list(),
[list(set(bond_df[bond_df["at1"]==x]["at2"].tolist() + \
bond_df[bond_df["at2"]==x]["at1"].tolist())) for x in df["id"]]))
template_name_after_reaction = "rxn1_ec_reacted.data_template"
build_react_template(df, bond_df, angle_df, dihedral_df, improper_df, selected_atoms, template_name_after_reaction)
viz_lmp_template(template_name_after_reaction)
Written to file: rxn1_ec_reacted.data_template
LAMMPS input script¶
T = 300 K, P = 1 atm, 27 Ethyl Cellulose Monomers going through
- 20,000 steps of Langevin+NPH with step size 0.5 fs;
- 10,000 steps of NVT with reaction enabled with step size 2.0 fs.
The reaction has a 100% happening when a pair of initiator atoms are detected with in 3.5 Å.
In [9]:
out_str = """variable Text equal 300.0 # reference temperature, in K
variable Ndump equal 1000
variable Ndump2 equal 100
units real
boundary p p p
atom_style full
include "opls.ff"
read_data "monomer.data" extra/special/per/atom 2
replicate 3 3 3
velocity all create ${Text} 4928459 dist gaussian
thermo ${Ndump}
thermo_style custom step dt cpu time temp press density pe ke xlo xhi ylo yhi zlo zhi atoms
thermo_modify flush yes
neigh_modify every 1 delay 0 check yes
variable padding equal 20
variable xLo equal xlo-${padding}
variable xHi equal xhi+${padding}
variable yLo equal ylo-${padding}
variable yHi equal yhi+${padding}
variable zLo equal zlo-${padding}
variable zHi equal zhi+${padding}
comm_style tiled
fix fxbalance all balance 100 1.1 rcb
# Adjust box dimensions:
change_box all x final ${xLo} ${xHi} units box
change_box all y final ${yLo} ${yHi} units box
change_box all z final ${zLo} ${zHi} units box
minimize 1.0e-20 1.0e-20 100000 100000
reset_timestep 0
# relax the cell
dump dump1 all custom ${Ndump} traj.lammpstrj.relax id mol type element mass x y z ix iy iz
dump_modify dump1 element """ + " ".join(mass_df["element"].to_list()) + """
timestep 2.0
fix fxlangevin all langevin ${Text} ${Text} $(150.0*dt) 12345
fix fxnve all nve
run 20000
write_data "monomers.pool.data"
unfix fxnve
unfix fxlangevin
undump dump1
# Reaction with NVT
timestep 2.0
dump dump2 all custom ${Ndump2} traj.lammpstrj id mol type element mass x y z ix iy iz
dump_modify dump2 element """ + " ".join(mass_df["element"].to_list()) + """
molecule mol1 """ + template_name_before_reaction + """
molecule mol2 """ + template_name_after_reaction + """
thermo ${Ndump2}
fix myrxns all bond/react stabilization yes statted_grp .03 react rxn1 all 1 0.0 3.5 mol1 mol2 """ + map_file_name + """
fix nvtReact statted_grp_REACT nvt temp ${Text} ${Text} $(100.0*dt)
fix velocityRescale bond_react_MASTER_group temp/rescale 1 ${Text} ${Text} 10 1
thermo_style custom step time temp press density f_myrxns[1] atoms
run 10000
write_data data.*
write_data data.final
"""
import io
with io.open("in.ec.stabilized", "w", newline="\n") as wf:
wf.write(out_str)
Now the simulation is ready to launch¶
In [10]:
!export OMP_NUM_THREADS=1 && ulimit -s unlimited && mpiexec -np 16 /home/xyan11/software/lmp20220623up4/build-cuda/lmp -in in.ec.stabilized
LAMMPS (23 Jun 2022 - Update 4)
using 1 OpenMP thread(s) per MPI task
Reading data file ...
orthogonal box = (-6.8 -1.2 -5.2) to (2.8 6.4 7.8)
2 by 2 by 4 MPI processor grid
reading atoms ...
36 atoms
scanning bonds ...
2 = max bonds/atom
scanning angles ...
6 = max angles/atom
scanning dihedrals ...
11 = max dihedrals/atom
scanning impropers ...
4 = max impropers/atom
reading bonds ...
36 bonds
reading angles ...
66 angles
reading dihedrals ...
90 dihedrals
reading impropers ...
20 impropers
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0.5
special bond factors coul: 0 0 0.5
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
16 = max # of 1-4 neighbors
21 = max # of special neighbors
special bonds CPU = 0.000 seconds
read_data CPU = 0.014 seconds
Replicating atoms ...
orthogonal box = (-6.8 -1.2 -5.2) to (22 21.6 33.8)
2 by 2 by 4 MPI processor grid
972 atoms
972 bonds
1782 angles
2430 dihedrals
540 impropers
Finding 1-2 1-3 1-4 neighbors ...
special bond factors lj: 0 0 0.5
special bond factors coul: 0 0 0.5
4 = max # of 1-2 neighbors
7 = max # of 1-3 neighbors
16 = max # of 1-4 neighbors
21 = max # of special neighbors
special bonds CPU = 0.001 seconds
replicate CPU = 0.008 seconds
Changing box ...
orthogonal box = (-26.8 -1.2 -5.2) to (42 21.6 33.8)
Changing box ...
orthogonal box = (-26.8 -21.2 -5.2) to (42 41.6 33.8)
Changing box ...
orthogonal box = (-26.8 -21.2 -25.2) to (42 41.6 53.8)
PPPM initialization ...
WARNING: System is not charge neutral, net charge = 0.0027 (src/kspace.cpp:327)
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.20258266
grid = 24 24 25
stencil order = 5
estimated absolute RMS force accuracy = 0.0028953053
estimated relative force accuracy = 8.719126e-06
using double precision FFTW3
3d grid and FFT values/proc = 3468 1152
Generated 820 of 820 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 10 9 12
1 neighbor lists, perpetual/occasional/extra = 1 0 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
Setting up cg style minimization ...
Unit style : real
Current step : 0
Per MPI rank memory allocation (min/avg/max) = 15.27 | 15.66 | 16.08 Mbytes
Step Dt CPU Time Temp Press Density PotEng KinEng Xlo Xhi Ylo Yhi Zlo Zhi Atoms
0 1 0 0 300 -8.8629696 0.031033818 633.17287 868.30997 -26.8 42 -21.2 41.6 -25.2 53.8 972
1000 1 1.8004668 1000 300 113.34444 0.031033818 -89.215268 868.30997 -26.8 42 -21.2 41.6 -25.2 53.8 972
1670 1 3.4328924 1670 300 112.94154 0.031033818 -136.34849 868.30997 -26.8 42 -21.2 41.6 -25.2 53.8 972
Loop time of 3.43327 on 16 procs for 1670 steps with 972 atoms
98.6% CPU use with 16 MPI tasks x 1 OpenMP threads
Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
633.172868877461 -136.34849172792 -136.34849172792
Force two-norm initial, final = 631.2321 4.5963283
Force max component initial, final = 45.21757 0.64485011
Final line search alpha, max atom move = 3.3702654e-10 2.173316e-10
Iterations, force evaluations = 1670 3390
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.00045239 | 0.65811 | 1.9084 | 89.3 | 19.17
Bond | 0.0009152 | 0.087443 | 0.2289 | 30.2 | 2.55
Kspace | 1.1428 | 2.5248 | 3.2858 | 51.5 | 73.54
Neigh | 0.019737 | 0.019961 | 0.020322 | 0.1 | 0.58
Comm | 0.065027 | 0.09075 | 0.111 | 6.3 | 2.64
Output | 4.9894e-05 | 7.0206e-05 | 0.00036154 | 0.0 | 0.00
Modify | 0 | 0 | 0 | 0.0 | 0.00
Other | | 0.05211 | | | 1.52
Nlocal: 60.75 ave 144 max 0 min
Histogram: 8 0 0 0 0 1 0 1 3 3
Nghost: 559.75 ave 855 max 287 min
Histogram: 8 0 0 0 0 0 0 0 1 7
Neighs: 14671.6 ave 42913 max 0 min
Histogram: 8 0 0 1 1 0 3 0 1 2
Total # of neighbors = 234745
Ave neighs/atom = 241.5072
Ave special neighs/atom = 10.5
Neighbor list builds = 37
Dangerous builds = 0
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.20258266
grid = 24 24 25
stencil order = 5
estimated absolute RMS force accuracy = 0.0028953053
estimated relative force accuracy = 8.719126e-06
using double precision FFTW3
3d grid and FFT values/proc = 3468 1152
Generated 820 of 820 mixed pair_coeff terms from geometric mixing rule
Setting up Verlet run ...
Unit style : real
Current step : 0
Time step : 2
Per MPI rank memory allocation (min/avg/max) = 15.58 | 16.01 | 16.48 Mbytes
Step Dt CPU Time Temp Press Density PotEng KinEng Xlo Xhi Ylo Yhi Zlo Zhi Atoms
0 2 0 0 300 112.94154 0.031033818 -136.34849 868.30997 -26.8 42 -21.2 41.6 -25.2 53.8 972
1000 2 0.6545845 2000 283.59939 65.514511 0.031033818 809.27501 820.84058 -26.8 42 -21.2 41.6 -25.2 53.8 972
2000 2 1.3078216 4000 300.24553 40.121457 0.031033818 829.68894 869.02061 -26.8 42 -21.2 41.6 -25.2 53.8 972
3000 2 1.9468943 6000 308.82944 -60.058771 0.031033818 823.75292 893.86559 -26.8 42 -21.2 41.6 -25.2 53.8 972
4000 2 2.6069884 8000 296.75892 -65.330246 0.031033818 799.46896 858.9291 -26.8 42 -21.2 41.6 -25.2 53.8 972
5000 2 3.252848 10000 303.58668 -17.699569 0.031033818 828.41948 878.69114 -26.8 42 -21.2 41.6 -25.2 53.8 972
6000 2 3.8597478 12000 297.36792 -1.6548949 0.031033818 841.5097 860.69178 -26.8 42 -21.2 41.6 -25.2 53.8 972
7000 2 4.5000472 14000 294.62181 -50.346715 0.031033818 787.49409 852.74351 -26.8 42 -21.2 41.6 -25.2 53.8 972
8000 2 5.142015 16000 301.71252 -112.91046 0.031033818 811.92558 873.26662 -26.8 42 -21.2 41.6 -25.2 53.8 972
9000 2 5.8237193 18000 308.94216 -1.656379 0.031033818 789.472 894.19185 -26.8 42 -21.2 41.6 -25.2 53.8 972
10000 2 6.4562026 20000 309.65983 76.398829 0.031033818 769.0346 896.26907 -26.8 42 -21.2 41.6 -25.2 53.8 972
11000 2 7.1089618 22000 311.60381 -87.261186 0.031033818 774.80169 901.89565 -26.8 42 -21.2 41.6 -25.2 53.8 972
12000 2 7.7740784 24000 301.65443 -12.548905 0.031033818 811.37676 873.09849 -26.8 42 -21.2 41.6 -25.2 53.8 972
13000 2 8.4463021 26000 316.39572 103.79187 0.031033818 824.66673 915.76519 -26.8 42 -21.2 41.6 -25.2 53.8 972
14000 2 9.0861241 28000 300.41137 46.270585 0.031033818 805.7395 869.50063 -26.8 42 -21.2 41.6 -25.2 53.8 972
15000 2 9.7524273 30000 288.02101 -193.4265 0.031033818 787.92829 833.63839 -26.8 42 -21.2 41.6 -25.2 53.8 972
16000 2 10.387209 32000 308.13764 -11.22038 0.031033818 798.33848 891.86329 -26.8 42 -21.2 41.6 -25.2 53.8 972
17000 2 11.072138 34000 293.43051 -68.678988 0.031033818 791.4712 849.29546 -26.8 42 -21.2 41.6 -25.2 53.8 972
18000 2 11.745767 36000 300.33916 1.8123468 0.031033818 754.22754 869.29163 -26.8 42 -21.2 41.6 -25.2 53.8 972
19000 2 12.406646 38000 304.51365 104.94022 0.031033818 790.16304 881.37412 -26.8 42 -21.2 41.6 -25.2 53.8 972
20000 2 13.044849 40000 289.22333 -184.2893 0.031033818 756.27092 837.11835 -26.8 42 -21.2 41.6 -25.2 53.8 972
Loop time of 13.045 on 16 procs for 20000 steps with 972 atoms
Performance: 264.929 ns/day, 0.091 hours/ns, 1533.154 timesteps/s
97.6% CPU use with 16 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.90437 | 2.836 | 3.9265 | 49.5 | 21.74
Bond | 0.45451 | 0.52822 | 0.59507 | 4.5 | 4.05
Kspace | 6.1771 | 7.4052 | 9.4653 | 33.8 | 56.77
Neigh | 0.42738 | 0.43926 | 0.4519 | 1.0 | 3.37
Comm | 1.1296 | 1.2874 | 1.506 | 11.7 | 9.87
Output | 0.022118 | 0.022326 | 0.025009 | 0.5 | 0.17
Modify | 0.39677 | 0.40034 | 0.40363 | 0.4 | 3.07
Other | | 0.1263 | | | 0.97
Nlocal: 60.75 ave 62 max 60 min
Histogram: 7 0 0 0 0 6 0 0 0 3
Nghost: 753.562 ave 878 max 564 min
Histogram: 1 3 3 1 0 0 0 0 0 8
Neighs: 14171.6 ave 19383 max 3662 min
Histogram: 1 1 1 0 1 1 2 3 1 5
Total # of neighbors = 226745
Ave neighs/atom = 233.27675
Ave special neighs/atom = 10.5
Neighbor list builds = 1379
Dangerous builds = 0
System init for write_data ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.20258266
grid = 24 24 25
stencil order = 5
estimated absolute RMS force accuracy = 0.0028953053
estimated relative force accuracy = 8.719126e-06
using double precision FFTW3
3d grid and FFT values/proc = 4860 1584
Generated 820 of 820 mixed pair_coeff terms from geometric mixing rule
Read molecule template mol1:
1 molecules
0 fragments
22 atoms with max type 29
20 bonds with max type 36
33 angles with max type 66
36 dihedrals with max type 89
10 impropers with max type 14
Read molecule template mol2:
1 molecules
0 fragments
22 atoms with max type 39
19 bonds with max type 39
32 angles with max type 70
36 dihedrals with max type 102
10 impropers with max type 14
dynamic group bond_react_MASTER_group defined
dynamic group statted_grp_REACT defined
WARNING: New thermo_style command, previous thermo_modify settings will be lost (src/output.cpp:904)
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
Your simulation uses code contributions which should be cited:
- fix bond/react: reacter.org
The log file lists these citations in BibTeX format.
CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.20258266
grid = 24 24 25
stencil order = 5
estimated absolute RMS force accuracy = 0.0028953053
estimated relative force accuracy = 8.719126e-06
using double precision FFTW3
3d grid and FFT values/proc = 4860 1584
Generated 820 of 820 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
update every 1 steps, delay 0 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 14
ghost atom cutoff = 14
binsize = 7, bins = 10 9 12
2 neighbor lists, perpetual/occasional/extra = 1 1 0
(1) pair lj/cut/coul/long, perpetual
attributes: half, newton on
pair build: half/bin/newton
stencil: half/bin/3d
bin: standard
(2) fix bond/react, occasional, copy from (1)
attributes: half, newton on
pair build: copy
stencil: none
bin: none
Setting up Verlet run ...
Unit style : real
Current step : 20000
Time step : 2
Per MPI rank memory allocation (min/avg/max) = 15.96 | 16.41 | 16.86 Mbytes
Step Time Temp Press Density f_myrxns[1] Atoms
20000 40000 289.22333 -184.2893 0.031033818 0 972
20100 40200 295.98625 -54.511648 0.031033818 0 972
20200 40400 298.11989 10.545549 0.030946176 1 969
20300 40600 291.41822 -25.694543 0.030946176 1 969
20400 40800 307.48353 27.106385 0.030946176 1 969
20500 41000 296.74622 36.713904 0.030946176 1 969
20600 41200 292.08604 23.039929 0.030946176 1 969
20700 41400 310.07704 50.450008 0.030858535 2 966
20800 41600 309.39629 -130.26994 0.030858535 2 966
20900 41800 304.45735 -131.40812 0.030770894 3 963
21000 42000 302.70032 -66.465065 0.030683253 4 960
21100 42200 296.66637 62.592419 0.030683253 4 960
21200 42400 297.85775 -26.219443 0.030683253 4 960
21300 42600 305.06705 -31.350749 0.030683253 4 960
21400 42800 298.47308 -36.173607 0.030683253 4 960
21500 43000 286.2572 11.286889 0.030683253 4 960
21600 43200 299.60746 -50.782826 0.030595612 5 957
21700 43400 301.79386 43.839688 0.030595612 5 957
21800 43600 305.87617 58.866287 0.030595612 5 957
21900 43800 286.00885 -14.04956 0.030595612 5 957
22000 44000 299.71152 77.613741 0.030507971 6 954
22100 44200 296.52791 -57.91458 0.030507971 6 954
22200 44400 302.19767 41.515013 0.030507971 6 954
22300 44600 286.2857 -140.49128 0.030507971 6 954
22400 44800 289.48453 -120.74076 0.030507971 6 954
22500 45000 306.38221 72.872366 0.030507971 6 954
22600 45200 320.54223 -2.6939923 0.030507971 6 954
22700 45400 306.50681 -19.373991 0.030420329 7 951
22800 45600 293.58002 24.673414 0.030420329 7 951
22900 45800 292.99665 3.3734624 0.030420329 7 951
23000 46000 283.68233 -68.938744 0.030420329 7 951
23100 46200 300.60638 -20.143097 0.030420329 7 951
23200 46400 305.3901 34.497133 0.030420329 7 951
23300 46600 298.65821 -17.700974 0.030420329 7 951
23400 46800 291.29385 -43.019826 0.030420329 7 951
23500 47000 291.53777 -78.903556 0.030420329 7 951
23600 47200 293.32735 45.926931 0.030420329 7 951
23700 47400 301.32633 -45.420213 0.030420329 7 951
23800 47600 308.39877 -63.611415 0.030420329 7 951
23900 47800 303.99618 -55.389737 0.030420329 7 951
24000 48000 295.14408 87.444387 0.030420329 7 951
24100 48200 296.89876 -93.246318 0.030420329 7 951
24200 48400 288.96838 -67.777401 0.030420329 7 951
24300 48600 290.93438 0.60378528 0.030420329 7 951
24400 48800 306.15046 31.100928 0.030420329 7 951
24500 49000 290.02075 18.492969 0.030420329 7 951
24600 49200 283.7607 12.25217 0.030420329 7 951
24700 49400 290.48281 -8.6403342 0.030420329 7 951
24800 49600 292.32201 38.592311 0.030420329 7 951
24900 49800 295.10302 -19.174346 0.030420329 7 951
25000 50000 276.53317 20.283103 0.030420329 7 951
25100 50200 286.42508 -16.230755 0.030420329 7 951
25200 50400 296.35416 -7.3621836 0.030420329 7 951
25300 50600 298.6115 -113.15535 0.030420329 7 951
25400 50800 298.07206 -45.843409 0.030332688 8 948
25500 51000 302.27055 -40.658546 0.030332688 8 948
25600 51200 292.42952 64.617125 0.030332688 8 948
25700 51400 281.63114 26.776393 0.030332688 8 948
25800 51600 292.48025 -48.31263 0.030332688 8 948
25900 51800 291.8906 62.084926 0.030332688 8 948
26000 52000 294.85678 14.111059 0.030332688 8 948
26100 52200 301.42051 3.1929241 0.030332688 8 948
26200 52400 284.06709 54.933763 0.030332688 8 948
26300 52600 298.84718 -109.47099 0.030332688 8 948
26400 52800 284.34178 -127.21862 0.030332688 8 948
26500 53000 296.36544 53.033968 0.030332688 8 948
26600 53200 288.51117 -41.492693 0.030332688 8 948
26700 53400 290.65741 -75.791795 0.030332688 8 948
26800 53600 282.33012 -64.641624 0.030332688 8 948
26900 53800 296.26217 18.437153 0.030332688 8 948
27000 54000 296.47698 -103.75234 0.030332688 8 948
27100 54200 297.23335 -9.7121337 0.030332688 8 948
27200 54400 295.45571 21.312435 0.030332688 8 948
27300 54600 285.88935 -4.0715574 0.030332688 8 948
27400 54800 283.58433 -23.789186 0.030332688 8 948
27500 55000 288.55117 83.325158 0.030332688 8 948
27600 55200 303.76952 1.0530389 0.030332688 8 948
27700 55400 295.38029 -27.102736 0.030332688 8 948
27800 55600 283.13427 -20.706221 0.030332688 8 948
27900 55800 279.01404 -38.774781 0.030332688 8 948
28000 56000 291.29063 -22.440435 0.030332688 8 948
28100 56200 294.09575 33.047792 0.030332688 8 948
28200 56400 303.6914 23.365713 0.030332688 8 948
28300 56600 306.83677 -9.1989978 0.030332688 8 948
28400 56800 289.83777 -52.926943 0.030332688 8 948
28500 57000 296.36805 -121.56544 0.030332688 8 948
28600 57200 292.28098 -32.695137 0.030332688 8 948
28700 57400 283.48553 32.646232 0.030332688 8 948
28800 57600 283.35026 -62.964664 0.030332688 8 948
28900 57800 293.84577 -148.69361 0.030332688 8 948
29000 58000 307.42555 24.305508 0.030332688 8 948
29100 58200 303.46263 -2.9106013 0.030332688 8 948
29200 58400 288.75175 -78.534012 0.030332688 8 948
29300 58600 290.62803 6.9585265 0.030332688 8 948
29400 58800 289.36775 -4.2108553 0.030332688 8 948
29500 59000 296.29667 10.613981 0.030332688 8 948
29600 59200 300.80576 84.835887 0.030332688 8 948
29700 59400 298.17714 40.17785 0.030332688 8 948
29800 59600 294.63556 -35.703098 0.030332688 8 948
29900 59800 292.68868 -35.800507 0.030332688 8 948
30000 60000 292.56691 -48.788968 0.030332688 8 948
Loop time of 8.33579 on 16 procs for 10000 steps with 948 atoms
Performance: 207.299 ns/day, 0.116 hours/ns, 1199.646 timesteps/s
97.6% CPU use with 16 MPI tasks x 1 OpenMP threads
MPI task timing breakdown:
Section | min time | avg time | max time |%varavg| %total
---------------------------------------------------------------
Pair | 0.42976 | 1.3716 | 2.1333 | 37.3 | 16.45
Bond | 0.23119 | 0.26679 | 0.31333 | 3.7 | 3.20
Kspace | 2.9116 | 3.8539 | 4.8481 | 24.8 | 46.23
Neigh | 0.22533 | 0.23458 | 0.24241 | 1.1 | 2.81
Comm | 0.5987 | 0.70635 | 0.88003 | 9.7 | 8.47
Output | 0.10687 | 0.10704 | 0.1093 | 0.2 | 1.28
Modify | 1.7447 | 1.7546 | 1.7631 | 0.4 | 21.05
Other | | 0.04093 | | | 0.49
Nlocal: 59.25 ave 66 max 52 min
Histogram: 2 1 2 1 1 2 1 3 0 3
Nghost: 866.562 ave 1044 max 640 min
Histogram: 2 0 0 2 2 6 0 0 0 4
Neighs: 13528.7 ave 22289 max 3211 min
Histogram: 2 1 0 0 3 4 1 1 3 1
Total # of neighbors = 216459
Ave neighs/atom = 228.33228
Ave special neighs/atom = 10.71519
Neighbor list builds = 741
Dangerous builds = 0
System init for write_data ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.20029103
grid = 24 24 24
stencil order = 5
estimated absolute RMS force accuracy = 0.0030385851
estimated relative force accuracy = 9.1506087e-06
using double precision FFTW3
3d grid and FFT values/proc = 6688 2240
Generated 820 of 820 mixed pair_coeff terms from geometric mixing rule
System init for write_data ...
PPPM initialization ...
using 12-bit tables for long-range coulomb (src/kspace.cpp:342)
G vector (1/distance) = 0.20029103
grid = 24 24 24
stencil order = 5
estimated absolute RMS force accuracy = 0.0030385851
estimated relative force accuracy = 9.1506087e-06
using double precision FFTW3
3d grid and FFT values/proc = 6688 2240
Generated 820 of 820 mixed pair_coeff terms from geometric mixing rule
Total wall time: 0:00:24
Analysis of the reaction¶
In [11]:
log_str = None
with io.open("log.lammps", "r", newline="\n") as rf:
log_str = rf.read()
import pandas as pd
import numpy as np
timestep = 1.0
react_df = pd.read_csv(io.StringIO(log_str.split("bytes")[-1].split("Loop time")[0]), sep=r"\s+", header=0)
react_df["react_instant"] = np.insert(np.diff(react_df["f_myrxns[1]"]), 0, 0)
react_df["time_ps"] = (react_df["Step"] - react_df.at[0, "Step"]) * timestep / 1000.
react_df
Out[11]:
| Step | Time | Temp | Press | Density | f_myrxns[1] | Atoms | react_instant | time_ps | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 20000 | 40000 | 289.22333 | -184.289300 | 0.031034 | 0 | 972 | 0 | 0.0 |
| 1 | 20100 | 40200 | 295.98625 | -54.511648 | 0.031034 | 0 | 972 | 0 | 0.1 |
| 2 | 20200 | 40400 | 298.11989 | 10.545549 | 0.030946 | 1 | 969 | 1 | 0.2 |
| 3 | 20300 | 40600 | 291.41822 | -25.694543 | 0.030946 | 1 | 969 | 0 | 0.3 |
| 4 | 20400 | 40800 | 307.48353 | 27.106385 | 0.030946 | 1 | 969 | 0 | 0.4 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 96 | 29600 | 59200 | 300.80576 | 84.835887 | 0.030333 | 8 | 948 | 0 | 9.6 |
| 97 | 29700 | 59400 | 298.17714 | 40.177850 | 0.030333 | 8 | 948 | 0 | 9.7 |
| 98 | 29800 | 59600 | 294.63556 | -35.703098 | 0.030333 | 8 | 948 | 0 | 9.8 |
| 99 | 29900 | 59800 | 292.68868 | -35.800507 | 0.030333 | 8 | 948 | 0 | 9.9 |
| 100 | 30000 | 60000 | 292.56691 | -48.788968 | 0.030333 | 8 | 948 | 0 | 10.0 |
101 rows × 9 columns
In [12]:
import plotly.graph_objects as go
from plotly.subplots import make_subplots
fig = make_subplots(specs=[[{"secondary_y": True}]])
fig.add_trace(
go.Scatter(x=react_df["time_ps"], y=react_df["react_instant"], name="Instantaneous Reaction Number"),
secondary_y=False,
)
fig.add_trace(
go.Scatter(x=react_df["time_ps"], y=react_df["f_myrxns[1]"], name="Cumulative Reaction Number"),
secondary_y=True,
)
fig.update_layout(
title_text="Number of Reactions vs. Time"
)
fig.update_xaxes(title_text="Reaction Time (ps)")
fig.update_yaxes(title_text="Instantaneous Reactions", secondary_y=False, showgrid=False)
fig.update_yaxes(title_text="Cumulative Reactions", secondary_y=True, showgrid=False)
fig.show()
The number of atoms in a reacted polymer molecule has a linear relation to the degree of polymerization (number of monomers)¶
$$N_{atoms}^{polymer} = (N_{dop}^{polymer} \times N_{atoms}^{monomer}) - (N_{dop}^{polymer} - 1) \times N_{atoms}^{water}$$ $$N_{atoms}^{polymer} = (N_{dop}^{polymer} \times 36) - (N_{dop}^{polymer} - 1) \times 3$$ $$N_{atoms}^{polymer} = N_{dop}^{polymer} \times 33 - 3$$ $$N_{dop}^{polymer} = \frac{N_{atoms}^{polymer} + 3}{36},\ \ \Big(N_{dop}^{polymer} \geq 2\ \ and\ \ N_{atoms}^{polymer} \geq 69\Big)$$ $$N_{dop}^{polymer} = 1,\ \ when\ \ N_{atoms}^{polymer} = 36$$ $$N_{dop}^{polymer} = \bigg\lfloor\frac{N_{atoms}^{polymer} + 3}{36}\bigg\rfloor$$
In [13]:
from polymaps.util import open_lmp_data, viz_mol
mass_df, df, bond_df, angle_df, dihedral_df, improper_df, \
[[xlo, xhi], [ylo, yhi], [zlo, zhi]], \
lj_coeff, bond_coeff, angle_coeff, dihedral_coeff, improper_coeff = open_lmp_data('data.final', viz=True, box=True, unwrap=True)
df
Out[13]:
| id | mol | type | q | x | y | z | ix | iy | iz | _x | _y | _z | element | color | size | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 40 | 2 | 4 | -0.6923 | -0.882380 | -3.804024 | 0.663180 | 0 | 0 | 0 | -0.882380 | -3.804024 | 0.663180 | O | rgba(255, 0, 0, 1.0) | 74.0 |
| 2 | 56 | 2 | 20 | 0.4110 | -0.970046 | -4.067368 | 1.599147 | 0 | 0 | 0 | -0.970046 | -4.067368 | 1.599147 | H | rgba(255, 255, 255, 1.0) | 46.0 |
| 3 | 60 | 2 | 24 | 0.0982 | -0.619815 | -5.260713 | 3.551134 | 0 | 0 | 0 | -0.619815 | -5.260713 | 3.551134 | H | rgba(255, 255, 255, 1.0) | 46.0 |
| 4 | 92 | 3 | 20 | 0.4110 | 0.682282 | -7.526284 | -0.698130 | 0 | 0 | 0 | 0.682282 | -7.526284 | -0.698130 | H | rgba(255, 255, 255, 1.0) | 46.0 |
| 5 | 53 | 2 | 17 | 0.4161 | -0.028561 | -4.153014 | -1.787572 | 0 | 0 | 0 | -0.028561 | -4.153014 | -1.787572 | H | rgba(255, 255, 255, 1.0) | 46.0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 944 | 815 | 18 | 23 | 0.0886 | 3.978399 | 0.663451 | 25.328454 | 0 | 0 | 0 | 3.978399 | 0.663451 | 25.328454 | H | rgba(255, 255, 255, 1.0) | 46.0 |
| 945 | 800 | 18 | 8 | -0.2477 | 3.697567 | -1.303694 | 26.077027 | 0 | 0 | 0 | 3.697567 | -1.303694 | 26.077027 | C | rgba(105, 105, 105, 1.0) | 77.0 |
| 946 | 836 | 19 | 8 | -0.2477 | 3.514176 | 7.115244 | 24.941209 | 0 | 0 | 0 | 3.514176 | 7.115244 | 24.941209 | C | rgba(105, 105, 105, 1.0) | 77.0 |
| 947 | 854 | 19 | 26 | 0.0982 | 3.454894 | 6.313397 | 25.618600 | 0 | 0 | 0 | 3.454894 | 6.313397 | 25.618600 | H | rgba(255, 255, 255, 1.0) | 46.0 |
| 948 | 853 | 19 | 25 | 0.0982 | 3.552899 | 7.972590 | 25.511356 | 0 | 0 | 0 | 3.552899 | 7.972590 | 25.511356 | H | rgba(255, 255, 255, 1.0) | 46.0 |
948 rows × 16 columns
In [15]:
from polymaps.util import map_discrete_colors
colormap_dict = map_discrete_colors(len(df["mol"].unique()))
df["color"] = df["mol"].map(colormap_dict)
viz_mol(df, bond_df)
In [ ]:
In [ ]: